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Abstract 

In this article we extend the exact simulation methods of Beskos et al. in to the solutions of 
one-dimensional stochastic differential equations involving the local time of the unknown process 
at point zero. In order to perform the method we compute the law of the skew Brownian motion 
with drift. The method presented in this article covers the case where the solution of the SDE 
with local time corresponds to a divergence form operator with a discontinuous coefficient at zero. 
Numerical examples are shown to illustrate the method and the performances are compared with 
more traditional discretization schemes. 
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1. Introduction 

1.1. Presentation 

Exact simulation methods for trajectories of one-dimensional SDEs has been a subject of much 
interest in the last years : see for example Q, 01) 29 1. Unlike the classical simulation 



methods which all involve some kind of discretization error (we mention 0| for the Euler Scheme) , the 
exact simulation methods are constructed in such a way that they do not present any discretization 
error under the strong hypothesis that the diffusion coefficient is constant and equal to one. In 
the last years, the original method presented in the fundamental article 0] has been extended to 
overcome various limitations of the initial algorithm ; it has been generalized to include the cases 
of unbounded drifts ([sl, j^). 

On another hand, the numerical simulation of SDEs corresponding to divergence form operators 

L = V. Qa(x)V 

involving a discontinuous coefficient a has been also the subject of various studies in the last years. 
Indeed these operators are of great importance since they appear in a wide range of modelling prob- 
lems involving diffusion phenomena in discontinuous media. Among applications, we can mention 



ecology (|9j|), geophysics (|l9|, |20|), astrophysics ([33|), or magneto/electroencephalography (tl3t]). 



In the one-dimensional context, various Random Walks and an Euler Scheme have been studied 



for the simulation of the solution of such SDEs : for Random Walks we mention ll] , [10| , |12l | , [1 



for the Euler Scheme see 22|, 23|, 2^ in the case where the discontinuity of the coefficient in the 
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divergence operator appears at point 0. Of course, for such SDEs, the order of discretization error 
of these discretization schemes is usually greater than those obtained in a more classical context. 

An important problem comes from the fact that SDEs corresponding to divergence form opera- 
tors do not enter the classical scope of SDEs covered by the exact simulation methods. 

The main difficulty is that these SDEs include an additional term, which involves in dimension 
one the local time of the unknown process (in dimension greater than one, it involves the local 
time of a one-dimensional auxiliary process; see In fact, the laws of the solution of such 

one-dimensional SDEs are no longer absolutely continuous with respect to the Wiener measure. 

In this paper we present a first attempt for the adaptation of the exact simulation methods of 
Q to one-dimensional SDEs with an additional term that involves the local time of the unknown 
process at point 0. Namely, our object of study is {Xt)t>o solution of 

dXt = a{Xt)dWt + biXt)dt + l3dL^t{X), Xo = x, (1) 

where < 1 (with /3 7^ 0) and Lf{X) is the symmetric local time of X in zero at time t ; the 
reason why we only deal with < 1 is that there is no solution to ([!]) when /3 > 1 (the case where 

= 1 corresponds to a refiected diffusion and we do not include it in our discussion). The reason 
why we restrict ourselves to the case /3 7^ is made clear in Section 3 (see Remark |2] and also the 
conclusion of the paper). 

Note that when a is identically equal to 1 and b is identically equal to 0, the solution (Xf) of ([T]) 
is a standard Skew Brownian Motion (SBM in short). Under mild assumptions concerning b and 
standard ellipticity conditions on a, it is known that there exists a unique strong solution {Xt)t>o 



to ([TJ (see [l6| for details). 



Let us emphasize that this work includes the situation where b ma y b e discontinuous at 0. So 



that the results of this paper are also suited for the situation stated in j23|, |2j| where the solution 
of ([1]) corresponds to a divergence form operator whose coefficient is discontinuous at (and is 
sufficiently smooth elsewhere). We show a numerical example to illustrate this interesting case. 

Let us now briefly explain our main idea. When o" = 1, we show that the law of {Xt)t>o (solution 
of ([T|)) is absolutely continuous with respect to the law of some Skew Brownian Motion (SBM) with 
a drift component. The reason why the SBM with drift appears naturally in our computations is 
explained in Section |3] (see Remark [T]) . 

So, contrary to the already mentioned discretization schemes where the standard SBM is used in 
force, we do not longer deal with a simple SBM but with a SBM that possesses a drift component. 
As a consequence, in order to adapt the method of [4] in this setting, we have to be able to simulate 
bridges of the SBM with drift. 

In the last section of the paper, we discuss the limitations of the initial algorithm. 

The main issue is to relax the boundedness assumptions made on the drift function b, as is done 
in jsj for "classical" SDEs. In jsl, the authors use some kind of factorisations for the sample state 
space of the standard Brownian Bridge, which are consequences of William's decomposition theorem 
for Brownian Motion. Proving similar factorisations for the Skew Brownian Bridge with drift seems 
difficult to us. Nevertheless, we have been able to apply a result stated in Pitman- Yor (2^ in the 
case of the standard Skew Brownian Bridge, which gives a first partial result. Unfortunately, we 
have not been able to relax the boundedness assumption on the drift function b and we think that 
much remains to do in this direction. 

In our concluding remarks we discuss the particular problem of being able to produce an exact 
simulation algorithm in the case /3 = 0, which in our opinion should be regarded as a special separate 
problem. We also draw bold lines for further investigation. 

1.2. Organisation of the paper 

In Section [2] we precise the hypotheses and define the problem we will deal with. We also 
introduce notations used in the sequel. In Section |3] we present the exact simulation algorithm. 
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adapted from Beskos et al. to our situation. It turns out that, in order to use the algorithm, 
we need to sample bridges of a Skew Brownian Motion with drift. Section U] is devoted to the 
computation of the transition probability density of the Skew Brownian Motion with drift. Then 
Section [5] explains how to sample bridges of a Skew Brownian Motion with drift, using rejection 
sampling with brownian bridges as proposals. Section [6] presents numerical experiments, including 
a divergence form case. Finally, we discuss possible extensions in Section [T] 



2. Exposition of the exact simulation problem. 

2.1. Exact simulation problem and first assumptions 

Denote C = C([0,T],M) the set of continuous mappings from [0,T] to R and C the Borel cr- field 
on C induced by the supreme norm. 

Let P be a probability measure on (C, C) and W a Brownian motion under P together with its 
completed natural filtration {J~t)t>o- We will denote = P (• | Wq = x). 

Throughout the whole paper, we will make the following assumptions 

- |/3|<1. 

— The function 6 : M — t- R is bounded and differentiable on R*'"*" and R*'^ with a possi- 
ble discontinuity at point {0}. We suppose that both limits lim^_s>o+ K-^) '■— ^(0+) ^^^^ 
limz^O-b{z) := b{0—) exist and are finite. The value 6(0) of the function 6 at is of no 
importance and can be fixed arbitrarily to some constant (possibly different from either 6(0+) 
or 6(0-)). 

We seek for an exact simulation algorithm of the paths of the solution of the one-dimensional 
Stochastic Differential Equation 

dXt = dWt + b{Xt)dt + pdL^tiX), Xo = x, (2) 

where L^{X) is the symmetric local time of X in zero at time t. 

2.2. Some recalls on SDEs of type ([2]) 
2.2.1. Existence and uniqueness 

Under the assumptions of the previous section, the equation ([2]) possesses a unique strong 
solution. In fact, performing the bijective change of variable g{x) i— t- (1 — j3)x\x>Q + (1 + (3)x\x<o 
allows to consider Yt := g{Xt) solution of a new transformed equation 

dYt = l/{gi^y{Ys)dWs + bog-^/{g^^y{Y,)ds 

without local time. Here (g^^)' denotes the half sum of the right and left derivatives of g^^. 

Since (g^^)' is bounded from below by a strictly positive constant, this transformed eq uation 



makes sense and well-known results for classical one dimensional SDEs (see for example |14l | Chap 
5. Section 5.5) ensure that under the assumptions of Subsection 12. H it possesses a weak solution. 
Then, g~^{Y) gives a weak solution of ([2]). 

The difficulty concerns strong uniqueness. Strong uniqueness for the solutions of equation ([2|) 
is proved with a direct application of Theorem 1.3 p. 5 5 in the fundamental article |lQ], which deals 
with a broader class of stochastic differential equations involving the local time of the unknown 
process. 

Note that when /3 = — lor/3 = +l, equation ([2|) possesses a unique strong solution, which is a 
reflected diffusion at 0, either reflected below (/3 = —1) or above (/3 = +1). 

Let us now briefly explain why there is no solution to equation ([2]) when > 1. Remember 
that we are working with the symmetric local time L^{X). Let us denote by L^'^{X) (resp. L^'\x)) 
the right-hand sided local time of the process X (resp. the left hand sided local time of X). It is an 
exercise to prove that if X is a solution of ([2]), then L^{X) = ^-^-L^'^ [X) and L^{X) = ^-^L^'\x). 
In particular, we see that when > 1 there is no solution to ([2]) (otherwise the symmetric local 
time of X would be negative !). 
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2.2.2. Strong Markov property 

The proof of the strong Markov property for solutions of equation ([2|) is a separate problem from 
the one of existence and uniqueness. We refer to jl5| for a rigorous proof. In the multidimensional 
context of diffusion processes with generalized drift, the proof of the Markov property may be found 
in 



32||. 



3. Exact simulation 

3.1. Notations and additional assumptions 
3.1.1. Notations 

Throughout the whole paper, we use the following notations : 



We note 



We set 



ao 2 



N\x) := / e~^dy. 

y/ ITT Jx 

M:=^S(0+)- 1^5(0-). (3) 

and define b{z) := h{z) — fi. 
We set 

z ^ (l){z) := y 

We set (j){z) := (p{z) — m with m = inf^giR 0(2;); the constant K denotes an upper bound of 
the function 6. 



We set B{u) := / b{y)dy, u € 
Jo 



^^.M will denote the SBM of parameter /3 and drift fi. That is to say B^'f^ is the strong 
solution of (|2]) in the case b = namely : 

dSf = dWt + ndt + (B^'i'). (4) 
We will denote p^'^^{t,x,y) the transition probability density of B^'^^. 

Note that, with this notation, p^'^^{t, x, y) is the transition probability density of the Brownian 
motion with constant drift /i G M, namely 

p 'f'{t,x,y) = -^exp{ 1. (5) 

Note also that p^'^{t,x,y) is the transition probability density of the SBM of parameter /3 
(without drift), see j2l| . 

We set 



tq := inf l^t > : B^'^ = Oj with the convention inf(0) = +oo. 
We will denote by h{x, .) the density of tq under P^. 
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3.1.2. Additional assumptions 

In this section, we will make the following additional assumptions : 

- /3 7^ (for the reason of this last hypothesis, see Remark [2|). 

- We assume that the function z i— )• 4'{z) is bounded. 

- We assume that the function u i— )• exp[S(u) — {u — x)^/2T] is integrable. 

3.2. Presentation of the exact simulation algorithm 
3.2.1. Change of probability 

Recall that in this case b{z) := b{z) — /i where /x is the constant defined by 

Note that since /3 7^ by assumption, this constant is well-defined. In the case where b is continuous 
at point {0}, observe that fi reduces to b{0). 
We have 

dXt = dWt + b{Xt)dt + ^dt + fidL^t{X). 
In particular, we may perform Girsanov's theorem (see Theorem 3.5.1 in fli'l) and we write 

dXt = dWf^ + fidt + /3dL° (X) , (6) 

where Wf^ := Wt + / b{Xs)ds is a Brownian motion under the new probability W"^^ defined by 
^0 

5^ = { £ "f^"'"'"''" - 1 [ '■'<^''*}' 

From our assumptions on 6, we are in position to apply the symmetric Ito-Tanaka formula to the 



function B{u) := / b{y)dy and {Xt)t>Q. 
Jo 

Applying the occupation's time formula, we obtain 



BiX.) - Bi.) ^ £ MM)±^.X. + i £ b'iX.n.^,,dt + M2±)_^,o 

(8) 



: r b{Xt)lx,^^dW^'' + r b{Xt)lx,^odt 
Jo Jo 



Jo \ 



where the last line comes from the definitions of b and (and the property dL^{X) = lxt=odL^{X)). 
From the fact that i{t G [0, T] : Xt = 0} = (where i stands for the Lebesgue measure), we see 
that ^ ^ ^ 

BiXr) - B{x) = [ b{Xt)dWf^ + fi [ b{Xt)dt + \ ( b'{Xt)dt. (9) 
Jo Jo ^ Jo 

Thus, ([7| implies that for any functional F{X) of the path up to time T, one has : 

Ep[F(X)] = E^sD [Fix) exp {^(Xt) - 5(x) - (l)iXt)dt}] , 

Jo 

where M.) = + ^""'^^ 
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Remark 1 Note that, because of the definition of b, there is no local time appearing in equality Q 
after the application of the Ito-Tanaka formula. This ensures that there is no local time involved 
in the exponential martingale of Girsanov's theorem, which makes it tractable for a numerical 
perspective. 

Retrospectively, this explains why in the sequel we have to deal with a Skew Brownian Motion 
with drift instead of a simple standard SBM. 

Remark 2 We now explain our assumption /3 7^ 0. 

Note that in the case /3 = 0, the constant fi is no more defined. In fact, in the case /3 = 0, 
because of the discontinuity of b, it is no longer possible to get rid of the local time as in (j8|). 
More precisely, there is no constant 6 such that proceeding as the computations in ([6]) and ([7|) with 
b{x) := b{x) — 9 we can cancel the local time term appearing in the exponential weight. 

For a more detailed discussion on the case /3 = the interested reader is invited to read the 
conclusion at the end of this paper. 

3.2.2. Exact simulation algorithm (after Beskos and al) 

Considering ([6|, we see that the law of X under W"^^ is given by p^'^{t,x,y)dy. 

Following the lines of Beskos et al. in |4,], and considering the computations performed in the 
above section, we give an algorithm that returns an exact simulation of a skeleton of (^t)ie[o,T] 
solution of (|2|) starting from xq. 



EXACT SIMULATION ALGORITHM FOR A SOLUTION OF 1^ starting from xq. 

1. Simulate a random variable Z according to the density 

h{y) = C exp {B{y) - i?(xo)) /''^(T, xq, y), 

WHERE C IS THE NORMALIZING CONSTANT SUCH THAT / h{y)dy = 1. KeEP IN MEMORY 
THE VALUE Z OF Z. 

2. Simulate a Poisson Point Process with unit density on [0,r] x [0,ii']. The result 

IS A random number n OF POINTS OF COORDINATES (ii, Zi), . . . , 2;„). 

3. Simulate a skeleton (sf^'^, . . . j-Sf^''^) conditioned on 5^'^ = xq and B^^'^ = z. 

4. If Vz g {1, . . . ,n} 4'{B^-^) < Zi accept the skeleton. Else return to step 1. 



This algorithm returns an exact sample of (Xt-i^, . . . , Xt„, Xt) (in particular we get an exact 
simulation of Xt, it is the value z oi Z used for an accepted trajectory). 

Note that in order to apply the methodology of j3| we have to be able to generate bridges of a 
drifted Skew Brownian Motion B^'^^. Indeed, this is the key one has to reach for in order to perform 
the Step 3. 

4. Computation of the law of the SBM with drift 

4.1. Recalls on the construction of the SBM using a "random flipping" of excursions 

In this paragraph, we present a construction of the Skew Brownian Motion - solution of (|4} 
with /X = and starting from x = - that gives an understanding of its relation with the standard 
Brownian motion. The construction is made out from a reflecting Brownian Motion with a change 
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of sign of each excursion with pro babihty (1 — /3) /2. It is explained in [28| page 487 exercise 2.16 
(we use the same notations as [28| in the explanations below). 

Suppose that x = 0. The construction is as follows : let (Yn)n>o be a sequence of independent 
r.v.'s taking the values 1 and —1 with probabilities (1 + /3) /2 and (1 — /3) /2 and independent of 
some Brownian Motion B. Let us set [tJ^^^^^ the natural filtration generated by B and satisfying 

the usual right continuous and completeness conditions and := Vt>o-^t^- 

We also denote Ti := cr{Yn : n > 0) the corresponding a-algebra generated by the whole sequence 
{Yn)n>o and £ := a{es '■ s > 0) the cr-algebra generated by all the excursions of B so that £ C J-^ . 
The "good" time clock for the point-wise excursion process (es)s>o is the process of time-change 
{jt)t>Q defined as the r.c.l inverse of the local time (L^ (i?))^^^ : so that the excursion process 
{^s)s>Q may be viewed as a Point Poisson Process on Co-5>o running in the local time scale. 

For each oj in the set on which B is defined, the set of excursions es(w) is countable and may be 
ordered. Define a process B^ by putting B^{ijj) = l^^^g) (a;)|es(t — Ts-{uj),u})\ if r^- ^ t < Tg and 
where is the ns(e)-th excursion in the above ordering. The random number ns{e) is a random 
variable measurable w.r.t. £, which depends on the whole excursion process e = {eu)u>o and the 
time variable s in the local time scale. It may be proved that the process thus obtained is a Markov 
process and that it is a Skew Brownian motion of parameter /3 starting at x = 0. 

By construction the sigma algebras Ti = cr(l^ : n > 0) and are independent ; in particular, 
if we denote by i?(e) the end point of excursion e, then we have Tt{uj) = R{^s{^)) and thus Tf is 

s<t 

measurable w.r.t. J^^ . Note that this construction implies that almost surely, L^{B^) = Lf{\B\) = 
L^{B) for any t > 0. So (^Lt{B^)^^^^ may be recovered as the r.c.l inverse of {Tt)t>o- Consequently, 
it is adapted to (•^t^)j>Q and it is independent of Ti. 

Remark 3 For a possible extension of this "random flipping of excursions" method for the con- 
struction of the solution to the more general equation ([2]) (and possibly solutions of ([1])), we refer 
to the article of Lejay fl7| . which gives a decomposition of the Ito measure associated to X in the 



general context of solutions of ([T]). At least in the context of equation ([2|), the result stated in |l7l | 
should allow to perform a construction along the same lines as above, flipping the excursions of 
some reflected process whose law should be the same as the law of \X\. However and up to our 
knowledge, such construction has never been explicitly written down in the literature, even in the 
context of equation ([2]). 

4-2. Computation of the joint law of SBM and its local time 

Let us begin with a direct consequence of the construction explained above in Subsection 14.11 

Lemma 1 Let W be a Brownian motion defined on (C, C,P) and B^ the strong solution of (j4| with 
/i = 0. 

We have for all t > 0, 

P°[|5f| e dy;L^t{B^) e dl] =P°[|Wf| e dy;L°(VF) e dl]. (10) 

Remark 4 We even have that the process {\B^\, Ll{Bl^))t>o is distributed as {\Wt\, L^{W))t>o 
under P°. This is a clear consequence of the construction explained in paragraph 14.11 Another 
way to prove this fact is to check that their common distribution is the one of {Mj^ — Wt, Mj^)t>Q 
under P*^, where Mj^ = maxo<s<t W^. This is related to the Levy theorem, as stated for instance 



in Theorem 3.6.17 in |l4| . where it is proved by using the Skorokhod method. 



Let us now state an intuitive result, which is somewhat not so easy to prove without using the 
construction explained in Subsection 14.11 The difficulty comes from the presence of the local time 
in the equalities below. 
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Lemma 2 We have for all t > 0, 



P°[Bf G dy;L',iB^) G dl] = ^P°[|i?f| G dy;L?(i?^) G d/] +:^P° [ -|sf | G rfy;L°(S^) G dl]. 

Proof. We start from the construction of SBM using a random flipping of excursions coming from 
a reflected Brownian Motion as explained at paragraph 14.11 

Let S be the space of real sequences {ak)kGn and denote <I> : 5 x N — t- M the coordinate function 
defined by $((afc)fcgN; 'i-) = On- From the independence of Ti and and since (L^(S'^))^^q is 
adapted w.r.t. (-^t^)j>Q, from the properties of the conditional expectation, 



l-/3™or 



^°[5f G dy;L°t{B'^) e dl] = ly>oE^ P°[sf G dy-L^t{B^) e dl\ 
lP°[yn,(e) >0;|e,(t-r,_H,a;)| Gdy;LO(i?'^) Gd/|^^] 
pO[$((n)fceNM,n,(eM)) >0| .F^ ] ; |e,(t - t,_H, w)| G dy; lO(S^) G 
^\Yn > 0] U=„4eM);|es(t-T,_(a;),w)| G dy-L%B^) G 



= ^F°[|i?f|Gdy;L°(i?'')Gd/]. 
Proceeding similarly on W_ we get, 



^°[sf G dy;L^t{B^) G d/" 



[-|sf|Gdy;L0(i?'5)Gd/]. 



therefore the result. 

Using the two last lemmas we can prove the following result. 



□ 



Proposition 1 Let W be a Brownian motion defined on (C, C,P) and B^ the strong solution of 
(jlD with /i = 0. 



We have for all t > 0, x > 0, 
P^[5f G dy-L^t{BP) G dl] 



+l,>o-^(exp{-(rf} - exp{-M£)!})dy5o(rf/) 



+lj^<olz>o- --i^ V 2t 



exp{ 



^dydl. 



Remark 5 The result of Proposition [T] appears as a corollary of a more general result stated in 
jsH- Note also that the result of Proposition [T] (and consequently the result stated in Proposition |2] 
in the next Section) differ slightly from results published by T. Appuhamillage et al. in the recent 
article [1] where there is a computational error (see also |2| for a discussion). We will detail the 
computations for the sake of completeness and clarification. 

Proof of Proposition [H Step 1. Combining the results of the Lemmas [T] and [2] we have 



Lj/>0^ 



^0[sf G d2/;LO(S/3) G d/] = i±^pO[|5f| G (iy;L0(S/3) G d/] 

i+/3tooJ|^^| g dy;LPi{W) e dl]. 
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Step 2. Let x > 0. As lj^>ol;>o P^[Sf G dy;L^{B^) e dl;t < tq] = 0, we have, using the strong 
Markov property, 

lj/>oli>oP^[sf e dy-Vl{BP) edl] = l2,>ol;>oP^[Bf G dy;L^t{Bl^) G cZ/;i > tq] 

= l^>oli>oE-[l{t>,„}P-[ i?f G dy; LO(i?/^) G dl \J^ro]] 

= iy>oii>o /o* P°[-Bf_3 e rfy; e dl ]hix, s)ds, 

where h{x, .) is the density of tq under P^. But h{x, .) is also the density of Tq = inf > : Wt = 0). 
And using the first step of the proof we have 

ly>ol/>oP°[5f_, G dy-Lts{B^) e dl] = l,>o^^P°[|m-.| e dy-Lts{W) G dZ]. 
Using again the strong Markov property wc get 
lj,>ol,>oP"[5f G dy;L^t{B^) G dl] = l,>oE-[l{,>To}i±^P^[ 1^1 € d2/;L?(W^) G |^To]] 

= li>o^^P"[|m| e dy;LO(W^) G > Tq] 
= li>o^^P"[|m|edy;LO(M/)Gdl] 

= lj,>ol;>o 7== exp I jdydl. 

Step 3. It is a consequence of the reflection principle that 

lj,>oP"[5f G dy-L'liBP) = 0] = lj,>oP^[i3f G (iy;^ < tq] 

1 / (y — x)^ (y + 
= l,,>n , — = I expi ) — expi — 



'^-'vwt — -^"p^ — ■ 

Using Step 1 to 3 we have the result for x > on R+ x M+. In order to retrieve the result on 
W_ X M+ we use Step 1 and 2 with replaced by and ly>o replaced by ly<Q, and the fact 
that for X > 0, ly<QW[Bl G dy-L^t{B^) = 0] = lj,<oP^[-Bf G < tq] = 0. □ 

-^.5. The law of the Skew Brownian motion with drift 
We have the following proposition. 
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Proposition 2 We have that, 



77= exp{/x(y - x) - \^iH} ( exp{ 



{y-x? 



2t + ^^{y -x)- 



2t 



}) 



X 



[1 - /3/xV2itexp{ ("+^+f^^ }Arc(M+|±^)] , 



i/x>0,y >0, 



1-/3 



2t 



+ fi{y-x) - 

[l — (3fi\/2TTt exp I 



2t 



I3fit+x—y 



if X > 0,y < 0. 



= exp{/i(y - x) - i^2i}(exp{- 



+ ^exp 



/2nt 



2t 



} - exp{- 



2t 



}) 



[l-/3//^/2^exp{ (-"-^+*^^)' 



2< 



i/x < 0,y < 0, 



2t 



, i/a;<0,y>0. 



+ ^x{y-x) - \^lH] 

[1 - /3/i^/2^exp { fa-^'+f }jv^( /3m^+|-^ )] ^ 



■ (N±M±1M!\ Arcr/9M+kl+|y|- 



Remark 6 It can be shown that the quantity 1 — l3fi\/2TTtexp — a }A^^(- ^ 

volved in p^'f^{t,x,y) remains strictly positive, whatever the sign of (see Remark [7]). 

Proof of Proposition [2l We have 

with Wl" = Wt + fit a Br ownian motion under defined by 



m- 



dF 



that under Q'^ the process B^''^ starting from is distributed as B^ starting from under 
For any bounded continuous function / and any t > 0, we have 



E|[/(Bf'^)] 



[/(Sf + x)] 
EO,[/(Bf + x) eM^^Wr - i/i^t}] 

/ 4^ fiy + x) expifiw - i/x2t}pO[5f g dy; Wt G dw] 
J J^, f{y) expifiw - yh}F^[B^ ^ dy; Wt - x G dw] 



(11) 



Suppose /3 > 0. 

We set ^xiz,l) = {z,z — x — /SI) which defines a bijection : M x 



Dr where 



{{y,w) £ R"^ : y - X > w}. Note that {B^,Wt - x) = ^^{B^ ,L^{B^)). Besides, almost surely, 
{B^, L^{Bl^)) G M X M+ and {B^ , W - x) £ D^. 

For X > 0, Proposition [T] ensures that the measure P^[sf e dy; L^{B^) G dl] has a density with 
respect to dydl on M x and gives mass to the segments of x {0} with the density 



;Bre<J!/;L?(B«) = 



/27rt 



2t 



(12) 



Let us denote A^; := {{y,w) £ x R : y = w + x} = $x(K+ x {0}). The measure P^[5f G 
dy; Wt — X £ dw] has a density yyiUj ^) with respect to dy dw on Dx \ = ^x(M ^ K*'"*")- But 
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it gives mass to the segments of the hne A^.. Let us denote ^^^{y,w) := {^i^{y,w), ^2^(11''^)) 
and notice that ^i^{y, w) = y. Let Ai C M+ and A = {{y, w) e M."^ : y e Ai, y = w + x} C A^.. As 
C M+ X {0} we have 

F-[{B^,Wt-x)(^A] = P-[(Bf,LO(i?^))G<I>-i(A)] 

= P-[i3f GcI>-i(^);L0(i?/3)=0] 

= P^'[sf G = 0]. 

Using this and (fT2]) in ([TT]) we get 

Ef[/(sf''^)] = Jf^^\^J-iy)eMf^w-'2l^H}F-[B^Gdy-Wt-xGdw] 
+ / /a. f(y) expifxw - lnH}F^[B^ £dy;Wt-x£ dw] 

= /r fiy) I-^ exp{fiw - \^iH}g'j^p^^{y, w)dw dy 

+ /(y)7bexp{^(y -rr) - - exp{-M£)!})dy. 

We now compute exp{/xti; — ^fi^t}g^fj ^)^'"^ with a change of variable and an integration 
by parts. We have for y > 0, 

fv-^ 1 „ (l + /3)(2^^f^ + a; + 2/) (li^+.+yf 
J exp{fiw--nt}g'^py^^{y,w)dw= ^ J e^""" -^== e dw. 



And, 



1^-- + ^ + j;)e- ' ' 2. ^' dw = /Je^^fe--) e-^^^-'K + x + y)e-^^^^^dw' 



which yields the desired result. The cases y < and /3 < are treated in a similar way. 

For the case x < 0, we perform the change of variable x — )• —x, y — )• —y, f3 — )• —f3 and /u — )• — 

□ 



5. Exact simulation of bridges of a Skew motion with drift 

For < t < T let us denote q^'^{t,T,a,h,y) the probability density of B^'^ knowing that 
Bq'^ = a and B^'^ = h. That is to say 

P[sf G dy I S^'^ = a , B^'f" = b] = q^'^^it, T, a, b, y)dy. 

Note that with these notations, q^'^(t,T,a,b,y) is the probability density of Wt knowing that 
Wq = a and Wt = b. As the law of the Brownian bridge is well known, sampling from q^'^{t, T, a, b, y) 
is easy. 

In order to sample along the law given by q^'f^{t,T,a,b,y) with a rejection algorithm using 
Brownian bridges values as proposals, we will use the two following results. 
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Lemma 3 Let a,b e M., < t < T . 
For G (—1,1) X M, we have 

VyGM, g^'^(t,T,a,6,y) = • (13) 

Proof. This comes from the Markov property for solutions of ([2]) (see Subsubsection l2".2.2l and the 
references therein). 

□ 

Lemma 4 We have for all € 

VO < t < T, Va, b,y£R, q°'^{t, T, a, b, y) = q^'°{t, T, a, b, y). 
Proof. It is a direct consequence of Equations ([5]) and (|13p . □ 
Lemma 5 Lei a, 6 e M, < t < T. 

- For {13, fi) G (—1, 1) X M with I3fi > 0, we have 

yy G R, q^'^it, T, a, b, y) < i^^'^^ q'^''>{t, T, a, b, y), (14) 

where 



and a := max(i^, ^-y^] 



For {/3,n) G (-1, 1) X M with < 0, set 



Y^''{t,z) := 1- /3fiV2Trtexp C 



Then, 
and 

In particular, 
where 



2t ' ' Vi 

p^'^{t,x,y) < 2a7^''^(i,|x|)/'^(i,x,y) (15) 

/•^(t,x,y) < 2aj^'^^{t,\y\)p^'>^it,x,y). (16) 

VyGM, q^'f'{t,T,a,b,y) <K^£^,q'''\t,T,a,b,y), (17) 



Proof. Case jSfi > 0. 

Let i > and x > 0. Looking at Proposition [2] it is clear that for y < 0, 

p^'^(t,x,y) < (l + /3)pO''^(t,x,y). 

For y > we have 

P^'''{t,x,y) < pO'^(t,x,y) + ^exp{-M! + ^(y-x)-i/i2i| 

< {l + ^3)p^^^^{t,x,y), 
12 



where we have used (y — x)^ < (y + x)^ (because x,y > 0). We can proceed in a similar way for 
X < and finally, we get that 

yt > 0, Vx, y G M, p^^^'it, X, y) < 2ap°'f'{t, x, y). (18) 

Thus, using the previous inequality gives 

pl^'>'{t,a,y)pP'f'{T -t,y,h) 



q^'^'it,T,a,b,y) 



pl^'^'{T,a,b) 
^ P^^'^jT, a, b) p°'^(t, g, y)pO-^(T - t, y, b) 
- " p/3.M(r,a,5) p<^'f'{T,a,b) ^ ' 

< o ' \i q'(t, T, a, b, y). 
pl^'i'{T,a,by ^ ' ' ' '^^ 



Case /3/x < 0. Let us denote T^'^'{t,x,y) := 1 - /j/x^/2i^ exp { (1^1+1^/^/^)' jNc^ ^f^t+\x\+\y\ y 
fixed X G M, y I— )• r^''^(t, x, y) is an even function. As we have 

2 /"OO 2 

Vz > 0, zeV / e~^du<l, (20) 

the function z i— ?■ \/27r exp(^)A^'^(z) has negative first derivative on M"*". Therefore y i— )• r^'^(i, x, y) 
is decreasing on M"*" and we have maXj^gK r^''^(i, x, y) = ^^'^{t, |x|). Using this and the same kind 
of computations than in the previous case we get (jlSp . As the roles of x and y are symmetric in 
T^'^{t,x,y), we get (116p . We then obtain (117p . using the same computations than for ([19}. □ 

Remark 7 Note that (I20p also allows to prove that p^'^{t,x,y) remains strictly positive (see Re- 
mark |6|) . 

Lemma |5]suggests to use the following rejection algorithm in order to sample along q^'^^{t, T, a, 6, y)dy 
(see for example Proposition 1 in [7]). 



EXACT SIMULATION ALGORITHM ALONG q^'''{t,T,a,b,y)dy 

1. Set K = K^'l^i^ if > 0, K = K^^^ ^^ otherwise. 

2. Sample y ALONG q°''^{t,Ta,b,y)dy. 

3. Evaluate 

^ 1 q^'^^{t,T,a,b,Y) 
>■ K q^fi{t,T,a,b,Y) " ' 

4. Simulate [/ ~ U{[Q,l\). If U < f{Y) accept the proposed value Y. Else return 
TO Step 2. 
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Exact Algorithm Bridges 



Acceptance Ratio 0.28 0.18 



Table 1: Acceptance ratios in Example 1. 



6. Numerical experiments 

Example 1. We first deal with a toy example. We consider the following SDE 

dXt = dWt-^ cos{^Xt)dt + /3dL^tiX), Xo = xo, (21) 

with /3 = 0.6, and xq = 0.2. Note that, here, the drift b{x) = — ^ cos(^a;) is bounded and of class 
C°° on the whole real line. 

For the exact procedure the constant drift involved in Subsection 13. 2| equals = 5(0) = — f- 
So we will have to sample bridges of SBM with non zero drift fi, using the results of Section [S] 

We have first to sample Xt from 

%) = Cexp{B{y)-B{xo))p^'''{T,xo,y) 

= Cexp(|(sin(fxo)-sin(fy)) -;u(y-:Eo))/'^(r,xo,y) 

(Step 1 of the Exact Algorithm). This can be done by rejecting standard normal random variables 
with mean xq. Indeed, using (|15p . we have here 



^ < 2a7^''^(r, |xo|)exp(5 - ^) /-"(T, xq, y). 

Then we accept or reject the proposed value Xt, using Steps 2 to 4 of the Exact Algorithm, with 
bridges of i?^''^, 

2 2 2 

(pix) = — cos (—X H sm(— j;) H , 

^ 8 ^5 ^ 20 ^5 ^ 20 

and K = ^g- as un upper bound for (f>. 

We plot on Figure [1] (top and bottom figures) the approximated density obtained with 10^ 
simulations of Xt, sampled with our exact procedure. On the top figure we plot the approximated 



densities obtained with 10 simulations of the Euler Scheme used in [22| and [23|, for decreasing 
time steps. We can observe the convergence of Euler type simulations to exact ones. Note that 
to have the Euler scheme fitting the exa ct p rocedure we have to take a fine time step (namely 



At = 10~^). This is because, as shown in 22|, the rate of weak convergence of the Euler scheme in 
this situation is of order (At)^/^~"^, for a smooth initial condition. 

On the bottom figure the approximated density is compared with the approximated densities 



obtained with 10^ simulations of the random walk based method studied in 12|, for decreasing 
space steps. Again we can observe the convergence of the process with discretization error. 

In Table [1] we report the empirical acceptance ratios for the rejection step using (p and the 
Poisson point process in the Exact Algorithm (this corresponds to the column Exact Algorithm 
in the table), and for the rejection sampling of bridges of the SBM with drift (this is the average 
acceptance ratio in this case). 

In Table |2] we report the CPU times needed to get the 10^ simulations, with the three different 
methods (and with the different discretization steps we have used). Programs were written in C- 
language and executed on a personal computer equipped with an Intel Core 2 duo processor, running 
at 2.23 Ghz. 
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Figure 1: Approximated densities of the positions at time T = 1.0 of 10^ paths of the solution of (|21ll starting from 
xo = 0.2: exact versus Euler with time step At = 10"", for n = 2,4 (top) and exact versus random walk with space 
steps h — j^, (bottom). 
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Exact 


Euler 

(At = 10"'', n = 2,4) 


Random Walk 
/ , 1 1 \ 
l/l = 10 5 200^ 


239s 


17s 
1680s 


3.52s 
1411s 



Table 2: CPU times for 10^ simulations of Xt- 



On this example the exact simulation is competitive, compared to schemes with very fine grids. 



Example 2. We want now to sample along the law of the continuous Markov process X 
generated by 

1 d d 

L = --ri'^-r ' I 

2 ax ax 



with 



a(x) 



(22) 



Note that a(0+) = 1 / 2 = a(O-). The coefficient a{x) is of class on M* " and M*'+, and 
uniformly strictly positive and bounded, which ensures the existence of X; in addition X solves 



(23) 



(see [18|, 11|). We define the Lamperti transformation ^(x) = dzj \J a(z) and set Yt ■= ^{Xt). 
Then 



^ V«(0+) + V a(O-) 



(24) 



(this follows from Proposition 3.1 in ll|; see also 18| and [25|). Firstly, note that 
1. Secondly, we have 

' 1 

{^)'{X) = 



V'a(0+)-y^(M 
Va(0+)+Va{0-) 



< 



if X > 



+ 6 ^y~fl^ ifx<0, 



2\/x2 + X + 1 - 2 



if X > 



and $ ^(y) = < 



-2-v/3x^ — X 



1 + 2^/2 if X < 0, 
As (y^)'(x) is bounded with bounded first derivative on 



^l+Vfa+2)2~3 



l-Vl-12[2-(v^-y/2)2] 



if y > 



if y < 0. 



and 



[>*,+ 



the explicitly known 



coefficients /3 = j ^^^^ ^(y) = |(\/^)' ° ""^(y) satisfy the assumptions of Subsub- 

section 13.1.21 Thus we can perform exact sampling from (|24p . and, applying the exact inverse 
transformation ^>~^, get samples from (|23p with absolutely no discretization error. 
Here we have, 

_ 1 o'(0+) - a'(O-) _ 26 
^ - 4y^(oT)-y^(Py " ~4(l-^/2)" 
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Figure 2: DIVERGENCE FORM OPERATOR: Approximated density of the positions at time T = 1 of lO'^ paths of 
the solution of H23p starting from a;o = 0.0: exact versus random walk with space step h = 3.10~^ and Euler scheme 
with At = 10"*. 



As we have 



-//y+ilog(^o$ !(?/)) 



if 2/ > 



B{y) 



-/xy + i[log(V^o ci>-i(y)) _ log(V2)] if y < 0, 



we can show (using again (|15p ). that for all yQ,y £ M and T > 0, 

exp(S(y) -i?(yo))/''^(T,yo,y) < y^y^2a7^''^(r, |yo|) e^^^V'^T, yo, 2/)- 

This allows to sample Yt from h{y) = Cexp{B{y) — B{yQ)) p^''^{T,yQ,y), by rejecting normal 
variables with mean yo and variance T. 

We then accept or reject the proposed value Yt by using bridges of B^''^ and 



((i/2)(VHy o ci.-i(y))^ + (i/2)((v^rv^) o ci>-i(y) 



with 



2x+l 



'4(a;2+x+l)3/2 (2x+l)Va;^+2:+l 
6a:- 1 3 



' (2x+l)i' 



if X > 



4(3x2 -a;+2)3/2 (2a;+l) V3x2-x+2 



+ 72 ^3x'-f;+^ if X < 0. 



We take = (^^~^/^) /4+(i4i-i/8)/2 upper bound for 0. We plot on Figure [2] the 

approximated density computed with 10^ simulations of Xt for xq = 0.0 and T = 1, obtained from 
the exact procedure. We plot on the same figure the approximated densities obtained with the 
Euler scheme and the random walk approximation mentioned in Example 1. 

We report in Table [3] the acceptance ratios. 
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Exact Algorithm Bridges 



Acceptance Ratio 



0.017 



0.5 



Table 3: Acceptance ratios in Example 2. 



Remark 8 Note that the acceptance ratio for the algorithm in the first example in Table [T] is quite 
low but decreases to less than 2% in Table [3] in the context of the second example. These figures 
are closely related to the measurement of the "distance" between the measure of the initial process 
from the reference measure and so these limitations of the algorithm arise even in the "classical" 
setting of the reference article ||4j| (for example with a rapidly varying drift). Nevertheless, in terms 
of CPU time, the performance of the algorithm seems quite competitive, in comparison with those 
of discretization schemes. 

Remark 9 Note that, at least graphically and contrary to what we can see on Figure [1] the 
transition density plotted on Figure |2] seems to be continuous at : this matches the well-known 
theoretical result, which asserts that the transition density of diffusion semigroups corresponding 
to elliptic divergence form operator of the form (|22p is always continuous. We refer to Stroock [30] 
for a proof based on the self-adjoint properties of these semi-groups and Nash's inequality. 

7. Discussion and concluding remarks 

7.1. An open problem : the path decomposition of a skew Brownian bridge 

An important issue for the extension of the initial exact simulation method is to overcome the 
restraining assumptions made on the drift function b (see Section [2.ip : namely, the assumption of 
boundedness for b. 

For example, it is frustrating that these assumptions do not allow us to simulate exactly what 
one may call the "Skewed Ornstein-Uhlenbeck" diffusion process. This difficulty appears even in 
the classical case (solutions of non skewed SDKs) and the fundamental reason is that we do not 
know how to simulate exactly a Poisson Point Process with c-finite intensity on the whole space 



In the classical case, where b is everywhere differentiable and no local time is involved (non 
skewed SDEs), this problem is solved by decomposing the trajectory of the standard Brownian 
bridge on [0, T] w.r.t. the space-time point where it attains its maximum or both its maximum and 
its minimum : we refer to [5j for a detailed presentation of this problem in the classical setting. 



Consequently, if one wants to overcome the restraining assumptions in Section 1271] concerning the 
drift function b, one has to search for such kind of decompositions for (at least) the Skew Brownian 
Motion (not to mention the drifted Skew Brownian Motion). Up to our knowledge, no results can 
be found in the literature concerning this decomposition and this open problem seems difficult to 
us. However, we give below some insight concerning this problem thanks to an application of a 
theoretical result stated in 1261 1. 



Let Tz'^ '■= inf(s > : Bs'^ = z). Set ux{x:,z) := ^'^^''^j which gives the Laplace 



transform of r^''^ at A > (with B^'^ = x). 

Proposition 3 (case fJ. = 0) 



Ml X M. 
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In the simple case where /i = 0, the function u\ is given by 
sinh (^V2X{z — x)^ 1 + /3 



ux{x;z) 



sinh(\/2Az) 2 cosh(^/2Az) - (1 - /3)e-^^ 
^ sinh(\/2Ax) 



sinh(V2Az) 

2\{x-z) 



2\x 



1 + /3 



2 cosh(V2Az) - (1 - /3)e-^2Az 



1-/3 



2cosh(V2Az) - (1 +/3)e 



2Az 



if z > X > 0, 

if X > z > 0, 
if X < < z, 

if > z > X, 
if z < < X, 



(25) 



sinh [ V 2X{z — x] 



1-/3 



smh(\/2Az) 2 cosh(\/2Az) - (1 + /3)e^^ 
^ sinh(\/2Ax) 



sinh(V2Az) 



if > X > z. 



Remark 10 Note that if /3 = 0, we retrieve after easy computations the well known result that 
gives the Laplace Transform of the law of the hitting time of z by a standard Brownian Motion 
starting from x. 

Proof. We only sketch the proof. The different cases may be easily conjectured from the description 
of the excursion measure for the SBM (Bs'^)s>o and the known facts concerning the standard 
Brownian Motion (decomposition of the different cases when a skew Brownian Motion reaches z 
starting from x). In order to check rigorously the validity of the result, one may verify that the 
formulas (|25p yield a solution of Dynkin's problem associated to the generator of {Bs'^)s>o namely : 



with 



1 (f 

2^"a(-;2;) = AnA(.;z) 
ux{z;z) = 1, 



nA(.;z) G{gGC0(IR)nC2((-oo,0)U(0,oo)) : (1 + /3)g'iO+) = (1 - f3)g'{0-)}. 



(26) 



□ 



A scale function s and the corresponding integrated speed measure m of a Skew Brownian 
Motion are given by 



■g^x if X > 
j^x if X < ' 



(/3 + l)x ifx>0 
ll-(3)x ifx<0. 



(see (21J). In particular, the density £^'^{t,x,y)dy of the SBM w.r.t. the speed measure m{dy) is 
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given by 



1 



exp{-<^}-..p{-^)) 



(1 + /3)V2^' 2t 
+_^exp{ -^^i^}, ifx>0,y>0; 



2t 



'2TTt 



1 ix — xi)^ 

exp { - — - — — }, ifx>0,y<Oorifx<0, y>0; 



'2TTt 



1 



2t 

e.p{-^}-»p{-^)) 



(1-/3)V2^' 2t 

1 :exp{-^^}, ifx<0,y<0 



2t 



(of course i^''^{t,x,y) = il^'°{t,y,x)). 
Let 

M^'° := sup i?f '0 ; p^/ := mf{s > : = Mf -O}. 

0<s<t 

Then, applying the results of Theorem 2 in Pitman- Yor (26| . we have the following proposition : 

Proposition 4 

1. For any a,b < z < oo, A > 0, we have that 



I \ p-'j 
■^{q-^Pt 1 



/3,0 




a, B, 



/3,0 
T 



uxia;z)uxiz;b) 



s{dz). 



2. Moreover, under . 



B, 



13,0 



a,B^/ 



^0 

'° : < s < n 
are independent, distributed respectively like 

^Sf '° : < s < rf linc/er P 



6, M^'*^ = z, /3y " = -u^ , i/ie path fragments 



(27) 



^T-s : < s < T - -u 



) 



5, 



/3,0 



a I given r. 



/3,0 



and 



B 



/3,0 



< s < r; 



/3,0 



under P 



/3,0 



/3,0 



An open problem is to find a description of these laws and to give a procedure in order to simulate 
these laws exactly. 

7.2. Concluding remarks 

In this paper we presented an extension of the exact simulation method of Q that permits 
to sample an exact skeleton of a one dimensional diffusion process skewed at 0. This method 
may be applied to diffusions related to strongly elliptic divergence form operators that possess a 
discontinuous coefficient at 0. The basic idea of this contribution depends highly on the possibility 
to perform a Girsanov transformation such that no local time appears in the Girsanov exponential 
weight and such that the reference measure is tractable. 

In our opinion, this first work should be extended in several directions. 

Firstly, it is necessary to give a complete treatment of the case /? = (see Remark [2]) . In this 
case, there still exists a way to perform a Girsanov transformation such that no local time appears 
in the Girsanov exponential weight^but then the reference measure becomes that of a Brownian 
motion with two-valued drift (see for an introduction to these particular types of Brownian 
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motions). As before, the difficulty arises for the simulation of the bridges in the Step 3 of the 
algorithm. 

Further digging shows that, in this particular situation, the solution of the exact simulation 
problem in the manner of jif] is closely related to the computation of joint laws for the position 
together with local and occupation times by an arbitrary Brownian bridge (with no drift) but 
conditioned on its final position and local time at 0. Even if there exists abundant litterature 
dealing with the Brownian bridge, there is no result for such joint laws. 

Consequently, we believe that performing a totally exact simulation algorithm in full generality 
in the case /3 = appears to be outside the scope of this paper. Note that a satisfactory treatment 
of the case /3 = is crucial if one has the objective to deal with the even more general case, where 
the discontinuity of b and the local time appear at distinct space points. 

Secondly, various questions arise in the treatment of "skewness" : how can we overcome the 
restraining boundedness assumption on the drift function b ? What about a one dimensional 
diffusion process skewed at a finite number of points ? 
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